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Abstract 

We discuss an Adams-type predictor-corrector method for the numer- 
ical solution of fractional differential equations. The method may be used 
both for linear and for nonlinear problems, and it may be extended to 
multi-term equations (involving more than one differential operator) too. 

1 Statement of the Problem 

In this paper we want to discuss an algorithm for the numerical solution of dif- 
ferential equations of fractional order, equipped with suitable initial conditions. 
To be precise, we first look at the fractional differential equation 

D?y{x) = f{x,y(x)), (1) 

where a > 0 (but not necessarily a £ N). We shall determine a solution of this 
equation on the interval [0,T], say, where T is a suitable positive number. At a 
later stage we will also investigate a more general problem (see §4) but for the 
moment we restrict ourselves to (1) in order to explain the fundamental ideas 
behind our strategy. 

It is well known [4, 31, 33] that there are many ways to define a fractional 
differential operator. The special operator £»? that we are using in (1) is defined 

D«y{x) = J m ~ a y {m) {x) (2) 
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where m := fa] is just the value a rounded up to the nearest integer, y ,rn) is 
the ordinary mth derivative of y (recall that m is an integer, so we are dealing 
with the classical situation), and 

J d z(x) = [ (x - tf-'ztfdt (3) 

1 \P) j 0 

is the Riemann-Liouville integral operator of order f3 > 0. It is common practice 
to call the operator D* the Caputo differential operator of order a because it 
has apparently first been used for the solution of practical problems by Caputo 
[5]. Note however that, according to [4, p. 11], it had already been introduced 
in a 19th century paper of Liouville. 

A typical feature of differential equations (classical or fractional) is that one 
needs to specify additional conditions to make sure that the solution is unique. 
In many situations these additional conditions describe certain properties of the 
solution at the beginning of the process, i.e. at the point x = 0. Therefore such 
a problem is called an initial value problem. It is easily seen that the number of 
initial conditions that one needs to specify in order to obtain a unique solution 
is m = \a \ . In particular, if 0 < a < 1 (which is the case in many applications), 
we have to specify just one condition. However the precise form of this condition 
is not arbitrary. Rather, when fractional differential equations are concerned, it 
turns out that there is a close connection between the type of the initial condition 
and the type of the fractional derivative. This is actually also the reason for us 
to choose the Caputo derivative and not the Riemann-Liouville derivative that is 
more commonly used in pure mathematics: For the Riemann-Liouville case, one 
would have to specify the values of certain fractional derivatives (and integrals) 
of the unknown solution at the initial point x = 0, cf. [33, §42]. However, when 
we are dealing with a concrete physical application then the unknown quantity 
y will typically have a certain physical meaning (e.g. a dislocation), but it is 
not clear what the physical meaning of a fractional derivative of y is, and hence 
it is also not clear how such a quantity can be measured. In other words, 
the required data simply will not be available in practice. When we deal with 
Caputo derivatives however, the situation is different. We may [11] specify the 
initial values 2/(0), y'(0), . . . , y (Tn ~ 1} (0), i.e. the function value itself and integer- 
order derivatives. These data typically have a well understood physical meaning 
and can be measured. 

Note that Lorenzo and Hartley [27] have discussed the problem of finding the 
correct form of the initial conditions in a more general setting (not necessarily 
assuming that the entire history of the process can be observed). 

Thus we combine our fractional differential equation (1) with initial condi- 
tions 

y (t) (0) = yj t k) , k = 0, 1, ■ ■ ■ ,m — 1, (4) 

where once again m = fa] and the real numbers y l 0 k) , k = 0, 1, . . . , m - 1, are 
assumed to be given. It then turns out that, under some very weak conditions 
on the function / on the right-hand side of the differential equation, we indeed 
can say that a solution exists and that this solution is uniquely determined [11]. 
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Specifically these conditions are (a) the continuity of / with respect to both its 
arguments and (b) a Lipschitz condition with respect to the second argument. 
We explicitly mention that these are the only conditions that we need to impose 
on /; in particular there is no need to assume that / is a linear function. As 
we shall see in §2, there is no need to impose a linearity assumption on / when 
we derive and discuss our algorithm. Thus we now have a tool that allows 
a convenient solution of nonlinear problems. Apparently the development of 
such nonlinear models has in the past been hampered by the lack of suitable 
(numerical or analytical) solution methods. 

We refer to [11] for a more thorough mathematical analysis of the initial 
value problem defined by eqs. (1) and (4). Some additional properties of the 
Caputo differential operator D“ are discussed in [22, 23, 31]. 

Many authors formally use Riemann-Liouville fractional derivatives, defined 

by m 

D°y{*) - 

where once again m = [c*], instead of Caputo derivatives. Typically those 
authors then require homogeneous initial conditions. It is known [31] that under 
those homogeneous conditions the equations with Riemann-Liouville operators 
are equivalent to those with Caputo operators. We chose the Caputo version 
because it allows us to specify inhomogeneous initial conditions too if this is 
desired. As we have described above, for the Riemann-Liouville approach this 
generalization is connected with major practical difficulties, cf. , e.g., [9, 14]. 

A few numerical methods for fractional differential equations have been pre- 
sented in the literature; cf., e.g., [2, 3, 6, 8, 15, 16, 17, 18, 19, 20, 22, 28, 29, 30, 
31, 32, 34, 37]. However many of these methods are essentially ad hoc methods 
for very specific types of differential equations (often just linear equations or 
even smaller classes). It is not clear if they can be generalized and how they 
behave when they are applied to an equation that does not fit into the class 
for which the algorithm has originally been derived. Our scheme, on the other 
hand, has been constructed and analyzed for the fully general set of equations 
without any special assumptions, and is easy to implement on a computer. We 
therefore believe that it will be of practical significance and helpful for solving 
a broad class of problems. 

2 The Predictor-Corrector Algorithm 

In this section we shall derive the fundamental algorithm that we have developed 
for the solution of initial value problems with Caputo derivatives. The algorithm 
is a generalization of the classical Adams-Bashforth-Moulton integrator that is 
well known for the numerical solution of first-order problems [24, 25]. 

Our approach is based on the analytical property that the initial value prob- 
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lem (1), (4) is equivalent to the Volterra integral equation 

vM = 'e + fxb) f <* - 

in the sense that a continuous function is a solution of the initial value problem 
if and only if it is a solution of (5). For a brief derivation of this equivalence we 
refer to [11, Lemma 2.3]. Note that the sum outside the integral on the right- 
hand side is completely determined by the initial values and hence is known. 

In order to motivate the construction of our numerical method, we shall first 
briefly recall the idea behind the classical one-step Adams-Bashforth-Moulton 
algorithm for first-order equations. So, for a start, we focus our attention on 
the well-known initial-value problem for the first-order differential equation 

Dy(x) = f(x,y(x)), ( 6a ) 

t/(0) = 3/0- ( 6b ) 

We assume the function / to be such that a unique solution exists on some 
interval [0,T], say. Following [24, §111.1], we suggest the use of the predictor- 
corrector technique of Adams where, for the sake of simplicity, we assume that 
we are working on a uniform grid {t n = nh : n = 0 , 1 ,..., AT} with some 
integer N and h := T/N. The basic idea is, assuming that we have already 
calculated approximations yh(tj) « y{tj) ( j = 1, 2 , . . . ,n), that we try to obtain 
the approximation by means of the equation 

y{t n +i) = y{t n ) + J /(*.!/(*))<**• ( 7 ) 

This equation follows upon integration of (6a) on the interval [t„,t n+ i]. Of 
course, we know neither of the expressions on the right-hand side of eq. (7) 
exactly, but we do have an approximation for y{t n ), namely yh{tn)> that we 
can use instead. The integral is then replaced by the two-point trapezoidal 
quadrature formula 

J g(z)dz « ( ff ( a ) + g(b )) , 

thus giving an equation for the unknown approximation y h (t n+ i), it being 

yh{t n +l) = 2 th(tn) + 2 f i^ n + l > V^n+1 ))) > 

where again we have to replace j j(t n ) and y(t n + 1 ) by their approximations y h {t n ) 
and y h {t n+ i), respectively. This yields the equation for the implicit one-step 
Adams-Moulton method, which is 

Vh{tn+ 1) = yh{tn) + ^ [fi^n 3 Vh (^n)) + /(^n+1 1 Vh Un+ 1 ))] * (^) 
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The problem with this equation is that the unknown quantity yu{tn+ 1 ) appears 
on both sides, and due to the nonlinear nature of the function /, we cannot solve 
for yh{tn+i) directly in general. Therefore, we may use eq. ( 8 ) in an iterative 
process, inserting a preliminary approximation for y h (t n +i) in the right-hand 
side in order to determine a better approximation that we can then use. 

The required preliminary approximation y^{t n -r 1 ), the so-called predictor, 
is obtained in a very similar way, only replacing the trapezoidal quadrature 
formula by the rectangle rule 

f g{z)dz « {b- a)g(a), 

J a 

giving the explicit (forward Euler or one-step Adams-Bashforth) method 

y^(^n+l) — Vhi^n) + hf{t n i yh{tn))' (^) 

It is well known [24, p. 372] that the process defined by eq. (9) and 


2 m(Wi) = Vh{tn) + 2 [fi^Vh^n)) + /(Wr^(Wl))] ’ 

known as the one-step Adams-Bashforth-Moulton technique, is convergent of 
order 2, i.e. 

max \y(t n ) ~ Vh{tn ) I = 0(h 2 ). 

n=l, 2 ,...,AT 

Moreover, this method behaves satisfactorily from the point of view of its nu- 
merical stability [25, Chap. IV]. It is said to be of the PECE (Predict, Evaluate, 
Correct, Evaluate) type because, in a concrete implementation, we would start 
by calculating the predictor in eq. (9), then evaluate /(< n +i > Vh (^n-fi))j use 
to calculate the corrector in eq. ( 10 ), and finally evaluate /(tn+i> 2 //i.Un+i))- 
This result is stored for future use in the next integration step. 

Having introduced this concept, we now try to carry over the essential ideas 
to the fractional-order problem with some unavoidable modifications. The key is 
to derive an equation similar to (7). Fortunately, such an equation is available, 
namely eq. (5). This equation looks somewhat different from eq. (7), because 
the range of integration now starts at 0 instead of t n . This is a consequence 
of the non-local structure of the fractional-order differential operators. This 
however does not cause major problems in our attempts to generalize the Adams 
method. We simply use the product trapezoidal quadrature formula to replace 
the integral, where nodes t 3 {j = 0 , 1 , . . . , n + 1 ) are taken with respect to the 
weight function (t n+x - •) a ~ 1 - In other words ’ we a PP^ the approximation 



z) a - l g{z)dz 


f 


(^n+1 


\a — 1 - 


9n+ 1 (*)d * 5 


( 11 ) 


where p n +i is the piecewise linear interpolant for g with nodes and knots chosen 
at the tj, j = 0 , 1 , 2 , . . . , n -f 1. Using standard techniques from quadrature 
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theory (cf. [13]), we find that we can write the integral on the right-hand side 
of (11) as 


tn + l 


(^n-rl 


- z)° 1 g n+l (z)dz 


h°_ 

a(a 4- 1) 


77,-)- 1 

^ ] Q'j y n±l9(tj) 

3— 0 


where 

f n a+1 — (n - a)(n 4 1) Q if J? = 0, 

a jn+ i= (n-i + 2)“ +1 + M)“ +1 -2(n-j + l) a+1 if 1 < j < n, (12) 

1^1 if j = n 4 1. 

This then gives us our corrector formula (i.e. the fractional variant of the one- 
step Adams-Moulton method), which is 


Vh 


Tal-l ,k 

^n+1 (k) 


h{tn+l) = 2^ “TT^° + 




fc=0 


4 


A! 

h a 


n»+ 2 ) s 


r(a + 2) 


/ Un+li 2//i (^n4l )) 


(13) 


^ a j t n+3 / (^j ^ Vh Uj))» 


where we have used the functional equation of the Gamma function twice (which 
yields r(cr)(a((a 4 1) — r(cK 4 2)) and the fact that Un-i-i 7 ™4i ^ ■ 

The remaining problem is the determination of the predictor formula that 

we require to calculate the value y£(t„+i). The idea we use to S eneralize the 
one-step Adams-Bashforth method is the same as the one described above for 
the Adams-Moulton technique: We replace the integral on the right-hand side 
of eq. (5) by the product rectangle rule 


(t 


n41 


1 g(z)dz % frj,rt41 g(tj 


where now 


bj, n +i = —{{n+l-j) a -{n- j) a ) 


(14) 


(see also [13]). Thus, the predicted value y£(t„+i) is determined by the frac- 
tional Adams-Bashforth method 

i«(w) = E* 4“ + E w<t.: »<*»• as) 


k = 0 


r(a) 


i— o 


Our basic algorithm, the fractional Adams-Bashforth-Moulton method, is com- 
pletely described now by eqs. (15) and (13) with the weights flyn+i anc i 
being defined according to (12) and (14), respectively. 

We have thus completed the description of our numerical algorithm. The 
mathematical analysis of this method in [12] shows that we may expect the error 

to behave as x 

max \y{tj) - Vh(tj)\ = 0(h p ) (16a) 
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where 


p — min(2, 1 + a) (16b) 

and the quantities h and N are related according to h = T/N> and T is the 
upper bound of the interval on which we are looking for the solution. The 
reason for this rather special form of the exponent p is that one can prove that 
p must be the minimum of the order of the corrector (which is 2 in our case) 
and the order of the predictor method (which is 1 here) plus the order of the 
differential operator (viz., a). This is a well known fact for PECE algorithms 
for first-order equations (the case a = 1). In view of the fact that 1 < p < ^ we 
have a satisfactory globally valid error bound. We note that, quite in contrast 
to the behaviour of the algorithm described in [8], the convergence order p of the 
Adams-Bashforth-Moulton scheme increases as a, the order of the differential 
equation, increases. 

As far as the stability properties of our algorithm are concerend, we note that 
it follows (unsing the general methods of Lubich [29]) that the stability prop- 
erties of this fractional Adams-Bashforth-Moulton scheme are at least as good 
as the corresponding properties of its counterpart for first-order equations, i.e. 
the classical second-order Adams-Bashforth-Moulton method. The properties 
of the latter are described in detail, e.g., in [25, Chap. IV]. 

For a more detailed investigation, including numerical examples, we refer to 
[12]. The remainder of this paper will be devoted to more practical aspects. 
Specifically in §3 we shall introduce some modifications that can enhance the 
efficiency of our scheme, and in §4 we extend the capabilities of our integrator 
so that it can handle a more general class of equations. Finally we present some 
numerical examples in §6, and for the convenience of the reader we also provide 
an appendix where the algorithm is stated once again, but in a pseudo-code 
type of notation. 


3 Modifications of the Algorithm 

We have completed the description of the basic algorithm and its most important 
properties. Now we want to address the question of whether it is possible to 
improve the performance of the algorithm in certain respects. There are indeed 
a few ways to do this, and we shall mention some here. It is easily seen that most 
of these modifications are independent of each other, so the user may implement 
almost any combination of them as required. 

3.1 Reduction of the arithmetic complexity 

First of all, there is a fundamental problem associated with all fractional differ- 
ential operators (not only the Caputo version that we look at here): In contrast 
to differential operators of integer order, fractional derivatives are not local op- 
erators. This means that we cannot calculate D°y{x), say, solely on the basis 
of function values of y taken from a neighbourhood of the point x where we 
work. Rather, we have to take the entire history of y (i.e. all function values 
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y(t) for 0 < t < x) into account. This is clearly exhibited in the definition of the 
operator, cf. eq. (2). As a matter of fact, this property is highly desirable from 
the physical point of view because it allows us to model phenomena with mem- 
ory effects. But when it comes to numerical work we find that this non-locality 
leads to a significantly higher computational effort: The arithmetic complexity 
of our algorithm with step size h is 0(/i~ 2 ), whereas a comparable algorithm for 
an integer-order initial value problem (thus involving a local operator) would 
only give rise to a 0(/i _1 ) complexity. A number of ways have been suggested 
to overcome this difficulty. The first idea seems to have been the fixed memory 
principle of Podlubny [31]. However a close mathematical analysis [13] reveals 
that the reduced complexity is achieved at the price of a significant loss in the 
order of accuracy of the method. A more promising approach seems to be the 
nested memory concept of Ford and Simpson [21] that may well be applied to 
our algorithm. This concept leads to an 0(h~ l log A -1 ) complexity (so almost 
all of the additional work introduced by the non-locality is removed), and it 
can be shown that this idea retains the order of accuracy of the underlying 
algorithm. For more details we refer to [21]. 

3.2 Additional corrector iterations 

Recall that in the case of a very stiff equation, we mentioned that the sta- 
bility properties of the Adams-Bashforth-Moulton integrator may not be suffi- 
cient. However, it is well known [24, 25] that the pure one-step Adams-Moulton 
method (i.e. the trapezoidal method) possesses extremely good stability proper- 
ties. These are spoiled in the Adams-Bashforth-Moulton approach only by the 
fact that in eq. (13) we cannot replace the predictor approximation on the 
right-hand side by the corrector value y n + 1 because, in general, we cannot solve 
that equation exactly any more. The idea is now to find a better approximation 
for the exact solution of that equation than the rather simple one obtained by 
applying just one functional corrector iteration with the predictor as a starting 
value. 

There are two main ways to achieve this goal. The first one is to use the 
value obtained by the first iteration (the first corrector step) as a new predictor 
and apply another corrector step. Obviously, this procedure can be iterated 
any given number of times, M say; the resulting method is called a P(EC) M E 
algorithm. In this way, we find a method that is “closer” to the pure Adams- 
Moulton technique, and therefore, its stability properties are also closer to the 
better properties of this method. Taking this idea to the extreme, we could 
even refrain from stopping after M iterations and continue to iterate until con- 
vergence is achieved. This would (theoretically) lead to even better stability, 
but the computational cost could be prohibitive, and it may even happen that 
(due to rounding errors) convergence could not be achieved numerically in finite 
precision. Moreover the iteration may not converge unless the step size h of 
the method is sufficiently small. It is possible to make the P(EC) M E approach 
more efficient by not using the same number M of corrector iterations in every 
step. In regions of higher stiffness, more steps may be taken in order to retain 
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stability and to keep the error under control; whereas, in regions where stiffness 
is not a problem, high accuracy may be achieved already with a small choice of 
M, thus speeding up the algorithm. 

In practice, one can implement this feature in such a way that the user can 
supply an upper bound for the number JVl of corrector steps to be taken. Setting 
this bound to 1 is then equivalent to switching off this modification completely. 
Moreover, the user can specify a tolerance s > 0 to the effect that the iteration 
is stopped if two consecutive steps give results that differ by less than this 
tolerance even if the maximum number of iterations has not been reached. 

In §4 we shall see that there also may exist other reasons for replacing the 
plain PECE structure by a more complex P(EC) M E version with some suitable 
M. 

As mentioned above, this is not the only possible approach to get a better 
approximation to the solution of the corrector equation. The second idea that 
we may use to enhance the stability of the method is to find another way of 
solving eq. (13) with y £ +1 replaced by 2/n+i- The most obvious way to do this 
would be to use a Newton iteration. As a starting value, we can still use the 
Adams-Bashforth predictor. Since it is known that Newton iteration converges 
faster (locally) than the simple functional iteration of the PECE process, we 
may expect to come very close to the pure Adams-Moulton method in just a 
few iterations. However, in order to implement Newton’s method, we need to 
work with the Jacobian of the right-hand side of the differential equation. This 
can also be a source of numerical problems, and may even lead to very long run 
times. We have, therefore, refrained from implementing and testing this so far, 
but a future extension in this way is possible. 

Note that, since we keep the Adams-Moulton formula as the basis of eq. (13), 
the convergence order of these modified algorithms remains unchanged. More- 
over, this modification of the algorithm does not alter the order of the arithmetic 
complexity in terms of the step size h. Only the stability is affected by these 
modifications. Further, we stress that this modification (unlike the P(EG) E 
scheme) does not avoid the problems that we shall encounter in certain of the 
situations mentioned in §4. 

3.3 Richardson extrapolation 

Whereas the suggestions above were constructed in order to improve the run 
time or the stability of the algorithm without changing the convergence be- 
haviour, we now turn our attention towards a possibile way of speeding up the 
convergence. As noted in [12], numerical evidence suggets that the error of our 
Adams scheme, taken to approximate the exact solution at a fixed point T > 0, 
possesses an asymptotic expansion of the form 

*<T) - »(T) = E c > hV + X>' 1 ' + “ < I7 > 

j= 1 J=1 
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where *i, k 2 and k z are certain constants depending only on the solution y and 
satisfying k z > max(2Aq , k 2 + O')- In practice it is unlikely that these constants 
can be determined explicitly, but as a rule of thumb one may often work with 
small or moderate values like h\ — 3 and k 2 = 5. We shall see below that 
knowledge of the value of k z is only necessary in order to derive an error bound 
but not in order to actually implement the convergence speedup procedure. 

Notice that the asymptotic expansion begins with an h“ term and continues 
with h 1+a for 1 < a < 3, whereas it begins with h 1+c \ followed by /i 2 , for 
0 < a < 1. 

Our belief in the truth of this conjecture is not only supported by numer- 
ical results but also by the theoretical analysis of de Hoog and Weiss [7, §5] 
who show that asymptotic expansions of this form hold if we use the fractional 
Adams-Moulton method (i.e. if we solve the corrector equation exactly) and that 
a similar expansion can be derived for the fractional Adams-Bashforth method 
(using the predictor as the final approximation rather than correcting once with 
the Adams-Moulton formula). Assuming that the conjecture is correct, we may 
use the relation (17) to eliminate the leading terms and thus obtain an approx- 
imation that converges faster than the original one. According to the usual 
constructions (cf., e.g., Walz [36]), this results in a Richardson extrapolation 
scheme that has the following form. We begin by sorting the exponents of h 
in eq. (17) in increasing order, thus obtaining a sequence j3> • • • where, in 

the case 0 < q < 1, we have j\ = 1 + <a, j 2 = 2, j z — 2 4- ol , j 4 = 3 4- ct, — 4, 
j 6 = 4 + a etc., whereas for 1 < a < 2 we have ji - 2, j 2 = 1 + a, j 3 = 2 + a, 
j 4 — 4 t j 5 = 3 + a, j 6 = 4 + a, etc. Obviously, we can proceed in a similar 
way for the case a > 2, but since this latter case seems to be of minor practical 
interest, we omit the details. Using this newly introduced terminology we can 
rewrite the asymptotic expansion (17) for the error in the form 


K - 1 

y(T)-y h (T)= 


k= 1 


( 18 ) 


with certain coefficients c£, some K € N and some positive real numbers j\ < 
j 2 < ... < j K% Note that, in a practical application, the j* are known by the 
considerations above whereas the coefficients c* k are in general unknown. 

The key observation is that we can use eq. (18) to improve the accuracy of our 
approximation even though we do not know those coefficients c k . We proceed 
by setting up a triangular array (a so-called Romberg tableau) of approximation 
values for y (T) of the form 


(0) 

Vh 



(0) 

Vh/2 

c 


(0) 
Vh / 4 

c 

(2) 
Vh/ 4 

(o) 

Vk/s 


ls> 

00 


(3) 

Vh/B 


(19) 
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The leftmost column of this array consists of the approximations determined by 
our Adams algorithm with the step sizes ft, ft/2, ft/4, etc.; i.e. = yn (T). We 
may terminate the succesive step size halving at any time we wish; this then 
determines the number of rows in the tableau. Having set up this (the zeroth) 
column, we turn our attention towards the next column. To determine these 
values, we use the asymptotic expansion and see that we may calculate a linear 
combination 


y { h % : = 




(o) 

Vh/ 2- 


in - 1 


Then, in view of (18), 


y(T) - 1 4/U 


oil (°) (0) 

2 3l y(T) — y(T) ^ ^hf 2 »~ l 

in - 1 2^ - 1 

2^ (y(T) - Vh,2*(T)) - (y(T) - y h ,»-i(T)) 

2TZ-1 ( E 4(2 + 0((2-»hy-<) j 
- (e 4(2^^ +o((2 i -^r)j 

2J1 — cT(2- |, /i) jl - - * — 

2jm — l 1V ' 2-?i — 1 

K - 1 

+ 53 4*(2-*‘/i) il + 0((2~^) JK ) 

fc=2 

K-l 

53 + 0((2-'*/i)«) 

k-2 


We can thus see that the ft jl term has been eliminated, and therefore the first 
column converges towards the true solution as O(hn) which is faster than the 
zeroth column. Moreover the error in the first column has got an asymptotic 
expansion that formally coincides with that of the zeroth column, except that 
the leading term is missing. Thus we may apply the same scheme again and 
again, and then we find that the i4h column is determined by the relation 


yh/2» 


V 


v yh/2» 


J-D 

yji/2*- 1 


in - 1 


There is almost no computational work involved in the calculation of the right 
part of the table: Only a few elementary operations need to be done. This 
is completely negligible when compared to the effort required for the zeroth 
column (where the basic Adams algorithm really needs to be carried out). It 
then turns out that the i/th column converges to y(T) with an error of the order 
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For the fractional differential equation algorithm presented in [8] (that is 
less suitable for nonlinear problems) a similar extrapolation strategy has been 
used successfully [16]. Note however that the exponents of the /i-terms in the 
asymptotic expansions for that algorithm are different from the exponents for 
our Adams scheme. 

We stress that, according to the results of Ford and Simpson [21], the extrap- 
olation procedure may also be applied (with the same exponents) if we replace 
our original Adams scheme by the modification mentioned in §3.1. 

3.4 General meshes 

Up to this point we have assumed the mesh points t 0 , 1 1 , . . . , t/v to be equispaced. 
In some practical applications it may be useful to relax this requirement. Indeed 
it is possible to do so. Of course we then find that the predictor weights b J<n+ 1 
and the corrector weights a ]n +i need to be changed. The corresponding expres- 
sions have been derived in [13, 14]. In this case the validity of the asymptotic 
expansion (and hence the possibility of working with Richardson extrapolation) 
will be destroyed, and in the error estimates the parameter h will then represent 
the maximum of the distances of consecutive mesh points. 


4 Multi- Term Fractional Differential Equations 

We now want to extend our algorithm to a more general class of equations. To 
this end we introduce the multi-term differential equation 


D° n y(x) = f(x,y(x),D^y{x),...,D°" l y(x)), (20a) 


equipped with initial conditions 

y ( fc )(0) = y { 0 k) , * = 0,l,...,m-l, (20b) 

where now m = [a„], and where we assume that the orders of the differential 
operators are sorted in such a way that 0 < or < < • * • < ot n . Since we 

have a total of n differential operators in (20a), we say that this is an n-term 
fractional differential equation. Obviously, we may recover the class of problems 
that we have considered so far upon setting n — 1 and a := a n . There are a 
number of instances where such multi-term problems arise; the earliest example 
seems to be the Bagley-Torvik equation 

D 2 .y{x) = bD 3 J 2 y{x) + q/(z) + f{x) 

that describes the motion of a rigid plate immersed into a Newtonian viscous 
fluid (cf. [35, 9]). Other special cases are Babenko’s model 

Dly{x) = F{x)D l , /2 y(x) + G{x)y(x) + H(x) 
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of a gas dissolved in a liquid [31, §8.3.3], or Koeller’s model 

D 2 m a y{x) = piD^y(x) + Pvy{x) + f{x) 

for a copolymer [26], where a is some fixed real number. 

Our strategy for the numerical solution of such initial value problems is 
based on the classical approach for higher-order differential equations: We want 
to rewrite the given problem as an equivalent system of equations, where each of 
those new equations contains only one differential operator and where all these 
differential operators are of the same order. It is easily oberved [10] however 
that this is possible only under some additional number-theoretic assumptions 
on the parameters ai, « 2 , • ■ ■ 5 &n- To be precise, we must assume them 

• either to be rational numbers 

• or to be commensurate real numbers (i.e. all the ratios cxj/ai - must be 
rational) contained in the interval (0, 1). 

Therefore we need to perform some preliminary manipulations before we can 
follow this path in the general case. Our algorithm thus consists of two stages. 

In the first stage we replace the given initial value problem (20) by the related 
problem 

D*'y(x) = f(x, y{x),D^ 1 y(x), , D°"~' y{x)), (21a) 

with initial conditions 

y( fc )(0) = yl k \ k - 0, - 1, (21b) 

with m = |’a n ). That is, we perturb the orders of the differential operators in 
the given equation, but all the other given data (the function / on the right-hand 
side and the initial values) remain unchanged. The perturbed orders di,...,a n 
are chosen according to the following criteria: 

(a) a, , . . . , Q n must be rational numbers, 

(b) r«„i — r^ni ) 

(c) gcd(l,di,. . . ,d„) should be as large as possible, 

(d) | cij - Qj\ should be as small as possible for all j. 

One can interpret this approach as a careful refinement of an idea sketched in 
[31, §9.2] in connection with applications in control theory. 

The essence of condition (a) is that it allows us to rewrite the newly con- 
structed initial value problem as the system of simple equations that we want 
to have. 

Condition (b) asserts that both the original problem (20) and the perturbed 
problem (21) require the same number of initial conditions. 
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In condition (c), we denote by gcd the greatest common divisor of its argu- 
ments. This concept is to be understood in the generalized sense here because 
the a 3 are not natural numbers. This means that 

gcd(zi, . . . ,z n ) = m&x{q E Q : z 3 /q E N V; = 1,2, . . .,n}, 

so that the greatest common divisor is the largest rational number that divides 
all the given numbers, and by “divides” we mean that the division gives a 
natural number as a result, without any remainder. As we shall see below, 
the size of this greatest common divisor essentially determines the size of the 
resulting system of equations and hence the computational complexity. 

Condition (d) assures y « y because, as shown in [10], 


I y - y lloo = o 


max i 

= 1,2,... ,n 


We can thus make sure that the solution y of the new problem (21) does not 
differ too much from the solution y of the original problem (20). 

It is obvious that the parameters a 1} . . . , d n are not determined uniquely by 
the four conditions, so there is some freedom for the user to trade off accuracy 
against speed. 

Note here that, to some extent, conditions (c) and (d) contradict each other: 
One can typically make the difference \a 3 - otj\ smaller by choosing a 3 as a 
rational number with a larger denominator, but this at the same time makes 
the greatest common divisor smaller. In a practical application one thus needs 
to balance these two conditions and find a useful compromise. 

In the second stage of our approach we want to solve the new initial value 
problem (21) numerically. In order to describe this second stage we introduce 
some notation. By 

7 := gcd(l,Qi,. . . ,5n) 

we denote the greatest common divisor (already encountered in condition (c)). 
Moreover we set 

A r := — . 

7 

Then we set up the system 


D?y 0 (cr) = 

Dly i(x) = j/2 (z)> 


DZy N _ 2 {x) 

DZvn_ a*) 

with initial conditions 


Vn-i(x), 

f(x, yo(x), y&i/-r( x )i ■ ■ • 


(°) = { 0 ° 


(22a) 
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for j 7 € No, 
else, 


(22b) 



and find [10] that this is equivalent to the perturbed initial value problem in 
the following sense: Assume that the solution of the system is the vector- valued 
function (y 0) . ■ -,yfr- J> and that y is the solution of (21). Then y 0 = y. 

We can immediately see that the dimension of the system is N, and by 
definition of this quantity we see that this is large if 7 is small and vice versa. 
This justifies condition (c) above. 

The initial conditions (22b) correspond to the given initial conditions in the 
sense that those components of (yo- ■ ■ ■ , yfi-i) that correspond to an integer- 
order derivative of y are assigned that initial value, whereas all other initial 
values are set to zero. A mathematical justification of this choice is given in [9]; 
a more physical argument can be found injl]. 

Then we note that we can rewrite this iV-dimensional problem in more com- 
pact vector notation as 

D:?(x) = F{x,Y{x)) (23a) 


with initial condition 

y{ o) - y 0 


(23b) 


where now V = • 'iVfj-. i)^ is & function of TV variables that maps to 


R 


1 N 


F(x,Y(x)) 


/ yi(x) \ 

V2{X) 


I Vn- i( x ) J 

'/(z>2/o(z).2/ai/7( x )> • • • > ya„- 1 /- r { x )) ' 


and Y 0 = (yo(0),yi(0), • . . ,2/^_ 1 (0)) T with the components of the vector being 
as in (22b). Formally the initial value problem (23) is identical with the problem 
consisting of eqs. (1) and (4) that we have considered in §§1, 2 and 3, the only 
difference being that we now consider vectors Y and F instead of scalars y 
and /. This difference produces no problems at all since we can simply apply 
the scalar scheme in a componentwise fashion to each component of the vector 
problem, and we find that all considerations made above (like the details of 
the implementation, the error analysis, and the possible modifications) remain 
valid. So we can now apply the Adams scheme to the vector problem and use 
the first component of the numerical solution Yf L say, viz. the approximation for 
the function y 0 (x), as an approximation for the desired solution y of the original 
multi-term fractional differential equation (20). 

Obviously, the error of this two-stage approximation scheme can be decom- 
posed additively into two parts, each of which can be attributed to one stage of 
the scheme, as 


y{x) - yo } h{ x ) = (y( x ) - y( x )) + (y( x ) ~ wA x )) = e ^( x ) + 

say, where ei is the error introduced in stage 1 by the perturbation of the given 
equation, and 62 is the error of stage 2 that is due to the fact that we can 
only solve the perturbed equation approximately and not exactly. Based on 
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this decomposition it may seem tempting to reduce the overall error simply 
by reducing the contribution from ei and leaving stage 2 unchanged. In view 
of what we said above, this can easily be achieved by choosing the atj values 
closer to the original ctj. However, we repeat that this is likely to increase the 
dimension of the system, and the very structure of the system (to be precise: the 
large number of zeros in the initial condition (22b)) has the unwanted effect of 
making this more precise system much more difficult to handle for a numerical 
integration scheme. Thus it will often be observed that the cavalier approach of 
improving in stage 1 and not changing stage 2 will lead to a small improvement 
in 6\ at the price of a major deterioration of €2 and hence a worse overall result. 


5 Further work 

One area that needs further investigation is this last dilemma. We have identified 
two possible ways forward, but the question remains open as to which of these 
(or some other approach) would be more appropriate. 

The first option we have identified is to decrease the step size h used in stage 
2 as the dimension N of the system increases (a useful rule of thumb would be 

to use h = c/N with some constant c). 

Our second approach would be to remove those zeros (that do not contribute 
to the overall solution, so that the numerical solution gets stuck at the initial 
value for the first approximately l/(2 7 ) steps) by replacing the plain PECE 
structure by a P(EC) M E version with a sufficiently large M (e.g. M - l/j). 

One could also combine these two ideas or use different values of M in 
different regions. For the present we can say that these two approaches both have 
the disadvantage of a significant increase in the computational effort, so one may 
often be better off with a quite crude approximation in the first stage followed 
by a simple version of the approximation algorithm in the second stage. This 
is particularly true if the given orders oi, . . . , a„ of the differential operators in 
(20a) are something like material constants that are known only up to a limited 
accuracy. 


6 Numerical Examples 


In this section we present some numerical examples to illustrate the error bounds 
derived above. We only considered examples where 0 < a < 2 since the case a > 
2 does not seem to be of major practical interest. In all cases, the computations 
were performed in the usual double precision arithmetic on a Pentium-based 
PC. 

Our first example deals with the equation 


D°y(x) 


J0320 8-« _ r(5 + q /2) 4— 0/2 + 9 r(Q + 1} 

r(9-a) T(5 — a/2) 4 

+ {^ x ° /2 ~ x 4 ) - [ y ( x )] 3/2 - 
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The initial conditions were chosen to be homogeneous (y(0) —0,y (0) 0, the 

latter only in the case a > 1). The exact solution of this initial value problem 

y(x)=x 8 -Zx 4+a ' 2 + \x a . 


We display some of the results in Tables 1 and 2. In particular, we have not only 
used our basic algorithm here, but we have also tried to increase the accuracy 
of the numerical solution by using the Richardson extrapolation idea of §3.3. In 
each case, the leftmost column shows the step size used; the following column 
gives the error of our scheme at a: = 1, and the columns after that give the 
extrapolated values. The bottom line (marked “EOC”) states the experimen- 
tally determined order of convergence for each of the columns on the right o 
the table. According to our theoretical considerations, these values should e 
1 + a 2, 2 + a, 3 + a, 4. 4 + a, ... in the case 0 < a < 1 and 2, 1 + a, 2 + a, 
4 3 + 4 + Qj . . . for 1 < a < 2. The numerical data in the following tables 

show that these values are reproduced approximately at least for a > 1 (see 
Table 1). In the case 0 < a < 1, displayed in Table 2, the situation seems to be 
less obvious. Apparently, we need to use much smaller values for h than m the 
case a > 1 before we can see that the asymptotic behaviour really sets in. t his 
would normally correspond to the situation that the coefficients of the leading 
terms are small in magnitude compared to the coefficients of the higher-or er 

terms > 

As usual, the notation -5.53(-3) stands for -5.53 • 10 3 , etc. 


Table 1: Errors for a = 1.25. 


step size 

error of 
Adams scheme 

extrapolated values 

1/10 

1/20 

1/40 

1/80 

1/160 

1/320 

1/640 

-5.53(-3) 
— 1.59(— 3) 
— 4.33( — 4) 
— 1.14( — 4) 
— 2.97(-5) 
— 7.66(-6) 
— 1.96(— 6) 

— 2.80(- 4) 
— 4.60(— 5) 
— 8.17(— 6) 
-1.54(-6) 
— 3.04(— 7) 
— 6.16(— 8) 

1.63( — 5) 
1.90(— 6) 
2.24(— 7) 
2.56(— 8) 
2.85(— 9) 

2.13(— 7) 
2.71(— 8) 
2.28(— 9) 
1.73(— 10) 

1.47(— 8) 
6.24( — 10) 
3.25(— 11) 

EOC 

1.97 

2.30 

3.17 

3.72 

4.26 


Our second example covers the case where the given function / (the right- 
hand side of the differential equation) is smooth. We look at the homogeneous 
linear equation 

D?y(i) = —y{x), y{ 0) = 1, 2/'(0) = 0 

(the second of the initial conditions only for a > 1 of course). It is well known 
that the exact solution is 

y(x) = E a (-X a ) 
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Table 2: Errors for a — 0.25. 


step size 

error of 
Adams scheme 


extrapolated 

values 


1/10 

1/20 

1/40 

1/80 

1/160 

1/320 

1/640 

2.50(-l) 
1 .81 ( — 2) 
3.61(-3) 
1.45(— 3) 
6.58(-4) 
2.97(-4) 
1.31(— 4) 

— 1.50( — 1) 
— 6.91 ( — 3) 
— 1.10( — 4) 
8.19(-5) 
3.49(— 5) 
1.12(— 5) 

4.09(— 2) 
2.16(— 3) - 

1.46(— 4) - 

1.92(— 5) - 

3.37(— 6) - 

8.15(— 3) 
3.89(-4) 
1.45( — 5) 
8.50(— 7) 

1.28(— 4) 
1.05(— 5) 
6.01( — 8) 

EOC 

1.18 

r 1.63 

2.51 

4.09 

7.44 

where 


oo 

_ k 




E °( Z ) ^r(afc + l) 

/c=0 v 


is the Mittag-Leffler function of order a. 

In Table 3 we state some numerical results for this problem in the case a < 1. 
The data given in the tables is the error of the Adams scheme at the point x = 1- 
We can see from the last line that the order of convergence is always close to 
1 + a. In contrast, Table 4 displays the case a > 1; here the results confirm the 
0{h 2 ) behaviour. This reflects the statement of eq. (16). 


Table 3: Errors for a < 1. 


h 

a = 0.1 

a = 0.3 

1/10 

-5.42(-3) 

— 1 .86( — 3) 

1/20 

— 1.22(-3) 

— 5.85( — 4) 

1/40 

— 4.40(— 4) 

— 1.97( — 4) 

1/80 

-1.68(-4) 

— 6.90(— 5) 

1/160 

-6.65(-5) 

— 2.49(— 5) 

1/320 

-2.68(-5) 

-9.18(-6) 

EOC 

1.31 

1.44 


a = 0.5 
— 1.30(— 3) 
— 3.93(— 4) 
— 1.26(— 4) 
— 4.18(— 5) 
— 1.42(— 5) 
— 4.86(— 6) 
1.54 


a = 0.7 

— 9.91(— 4) 
— 2.81(— 4) 
— 8.28(— 5) 
-2.50(— 5) 
-7.63(— 6) 

— 2.35(— 6) 

1.70 


a = 0.9 
— 7.51(— 4) 
— 1.91(— 4) 
— 4.99(— 5) 
— 1.32( — 5) 
— 3.54(— 6) 
— 9.48( — 7) 
1.90 


Appendix: Pseudo code of the algorithm 

Following the derivation and description of the Adams-Bashforth-Moulton al- 
gorithm in mathematical terms given in §2 we now state this algorithm in a 
pseudo code type notation. In this way the reader interested in implementing 
the method can easily do so in the language of his or her preference. 
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Table 4: Errors for a > 1. 


h 

a = 1.25 

a = 1.5 

a = 1.85 

1/10 

— 5.61(— 4) 

— 5.46(— 4) 

— 4.40(— 4) 

1/20 

-1.27(— 4) 

— 1.28(— 4) 

— 1.07( — 4) 

1/40 

— 2.90(— 5) 

— 3.04(— 5) 

— 2.65(— 5) 

1/80 

-6.68(-6) 

— 7.33(— 6) 

— 6.57(— 6) 

1/160 

— 1.55(— 6) 

— 1.78(— 6) 

— 1.63(— 6) 

1/320 

-3.63(-7) 

-4.37(— 7) 

— 4.07(— 7) 

EOC 

f 2.09 

2.03 

2.00 


Input variables 

/ - the real-valued function of two real variables that defines the right-hand 
side of the differential equation (1), 

a - the order of the differential equation (a positive real number), 

y 0 - an array of fa] real numbers that contains the initial values y(0), y'(0), 
J/ (r “ 1_1) (0), 

T - the upper bound of the interval where the solution is to be approximated 
(a positive real number), 

N - the number of time steps that the algorithm is to take (a positive integer). 

Output variables 

y - an array of N + 1 real numbers that contains the approximate solutions 
VT/nUT/N), j = 0, 1, . • • , N. We have ijt/nU t / n ) ~ v{ T / N ) where y is 
the exact solution of the given fractional differential equation. 


Internal variables 

h — the step size of the algorithm (a positive real number), 
m — the number of initial conditions specified (a positive integer), 
j, k — integer variables used as indices, 

a, f> - arrays of N + 1 real numbers that contain the weights of the corrector and 
predictor formulae, respectively, 

p - the predicted value (a real variable). 
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Comment 

In order to save memory, the arrays a and b are one-dimensional and not two- 
dimensional as indicated by the notation in §2. This is possible because the 
values a, k and bj k defined in eqs. (12) and (14), respectively, essentially have 
a convolution structure, i.e. they only depend on the difference k - j of the two 
indices. 


Body of the procedure 

h := T/N 
m := [a] 

FOR k := 1 TO N DO BEGIN 
6[Jfe] := k a — (k — l) a 
a[k] := ( k + 1)° +1 - 2 k a+1 + (k - 1) Q+1 
END 


y{0]:= 3 /o [0] 

FOR j := 1 TO N DO BEGIN 


771 — 1 


1 


p:- E + 


k~ 0 

771 — 1 


r(a + 1) t 


0 


*i := E ^voW 


Jc=0 


k\ 


T(a + 2) 


+ (O' - i) a+1 - 0 - 1 - «b' a ) /(°*»[°]) 


i - 1 


+ ^a[j - k]f(kh,y\k]) 


k—l 


END 
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